# uncomment following line to make an html
!jupyter nbconvert --to html CAESAR_Validation.ipynb

[NbConvertApp] Converting notebook CAESAR_Validation.ipynb to html
[NbConvertApp] Writing 7677883 bytes to CAESAR_Validation.html

import copy
import importlib
import caesar_utils
import math
import numpy as np
import netCDF4
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm, Normalize, ListedColormap
from metpy.units import pandas_dataframe_to_unit_arrays, units
from datetime import datetime, timedelta
from IPython.display import display
from bokeh.io import push_notebook, show, output_notebook
from bokeh.layouts import row
from bokeh.layouts import gridplot
from bokeh.plotting import figure, show
from bokeh.models import Title, CustomJS, Select, TextInput, Button, LinearAxis, Range1d
from bokeh.models.formatters import DatetimeTickFormatter
from bokeh.models.tickers import DatetimeTicker
from bokeh.palettes import Category10
import warnings
import itertools 
import holoviews as hv 
from holoviews import dim, opts
import hvplot.pandas
hv.extension('bokeh', 'matplotlib')
warnings.filterwarnings('ignore')
output_notebook()
%matplotlib inline

Loading BokehJS ...

#Maneuvers Found:
spd_run_times = {"tf01": [[69441, 70195]],
                 "tf02": [[73493, 73783]], #ncplot "utc seconds"...: "tf02": [[48293, 48583]],
                 "tf05": [[35378, 35914],  # directly over Scandinavian mountain range... 6000 ft
                          [42345, 42645]], # lee of mountains, 4000 ft
                 "rf01": [[46183, 46950]], #"rf01": [[20983, 21750]], # Speed over ocean
                 "rf04": [[37523, 38000]], #"rf04": [[12323, 15034]], # decel maneuver after descent, straight and level, further descent maneuver in turb. ocean
                 "rf06": [[39491, 39750]], # affected by icing, to be removed
                 "rf07": [[36671, 38159]], # accel and decel separate. Over mountains. 20k ft
                }
spd_change_times = {"tf02": [[58095, 58351]],
                    "tf04": [[47070, 47835]], # stepped speed change in Sweden with short turns.
                    "rf02": [[33348, 33698]], # accel from 115 to 133 m/s.
                    "rf03": [[9948, 10681]],  # stepped speed change with short turns. Not so level. Over the mountains.
                   }
yaw_times = {'tf01': [[70496,70590]], 
             'tf02': [[74035, 74139], 
                      [74205, 74342]], 
             'tf05': [[37739, 37804],  # accidental climb trimmed out, so short maneuver
                      [37936, 38033],  # good maneuver at 18k feet
                      [39951, 40034]]} # good maneuver at 10k feet
# pitch_times =   {"tf01": [[45444, 45556]],
#                  "tf02": [[49381, 49524],
#                           [49626, 49762]],
#                  "tf05": [[38240, 38385], # good pitch maneuver, 18k feet
#                           [40040, 40163]] # good pitch maneuver, 10k feet
#                 }
pitch_times =   {"tf01": [[70644, 70756]],
                 "tf02": [[74581, 74724],
                          [74826, 74962]],
                 "tf05": [[38240, 38385], # good pitch maneuver, 18k feet
                          [40040, 40163]] # good pitch maneuver, 10k feet
                }
reverse_hdg_times = {"tf01": [[45631, 45985], # heading 262
                              [46195, 46530]]} # heading 082
revrs_hdg_spd_change = {"tf04": [[47613, 47836], # reverse hdg maneuver, but with speed changes 136 vs 124 m/s
                                 [48027, 48349]]
                       }
racetrack_times = {"tf02": [[50010, 50186], # heading 311
                            [50294, 50411]], # heading 131
                   "tf05": [[39381, 39499],
                            [39606, 39705]]}
circle_times = {"tf01": [[46564, 46947], # 2+ left circles
                         [47003, 47348]], # 2+ right circles
                "tf02": [[51062, 51417], # 2+ left circles
                         [51445, 51807]], # 2+ right circles
                "tf04": [[48830, 49201], # 2+ right circles
                         [49221, 49665]], # ~3 left circles
                "tf06": [[46499, 46928], # 3- left circles
                         [47023, 47364]] # 2+ right circles (1 trimmed out due to brief unbank)
               }

def summarize(time_dict: dict) -> None:
    key_num = {k: np.array(v).shape[0] for k, v in time_dict.items()}
    total_maneuvers = np.sum([np.array(v).shape[0] for v in time_dict.values()])
    print(f"        Found in following flights: {key_num}")
    print(f"        Total maneuvers: {total_maneuvers}")

print("CAESAR Maneuver Summary:")
print("    Speed Runs:")
summarize(spd_run_times)

print("    Speed Change Maneuvers:")
summarize(spd_change_times)

print("    Yaw Maneuvers:")
summarize(yaw_times)

print("    Pitch Maneuvers:")
summarize(pitch_times)

print("    Reverse Heading Maneuvers:")
summarize(reverse_hdg_times)

print("    Circle Times:")
summarize(circle_times)

CAESAR Maneuver Summary:
    Speed Runs:
        Found in following flights: {'tf01': 1, 'tf02': 1, 'tf05': 2, 'rf01': 1, 'rf04': 1, 'rf06': 1, 'rf07': 1}
        Total maneuvers: 8
    Speed Change Maneuvers:
        Found in following flights: {'tf02': 1, 'tf04': 1, 'rf02': 1, 'rf03': 1}
        Total maneuvers: 4
    Yaw Maneuvers:
        Found in following flights: {'tf01': 1, 'tf02': 2, 'tf05': 3}
        Total maneuvers: 6
    Pitch Maneuvers:
        Found in following flights: {'tf01': 1, 'tf02': 2, 'tf05': 2}
        Total maneuvers: 5
    Reverse Heading Maneuvers:
        Found in following flights: {'tf01': 2}
        Total maneuvers: 2
    Circle Times:
        Found in following flights: {'tf01': 2, 'tf02': 2, 'tf04': 2, 'tf06': 2}
        Total maneuvers: 8

# Open all netcdf files and read them into data frames
importlib.reload(caesar_utils)

yr = '2024'
data_dir = '/Users/ckruse/qaqc/CAESAR/data/raf_data'
project = 'CAESAR'

# open all NetCDF files
nc_dict = caesar_utils.open_nc(data_dir) # Dictionary of flight NetCDFs

# read the variables selected in caesar_utils.py from each NetCDF file
data = {} # dictionary of DataFrame's
for flight in nc_dict.keys():
    data[flight] = caesar_utils.read_nc(nc_dict[flight])
    print(f"Done reading {flight}")

Found 6 ferry flights, 6 test flights, and 10 research flights
Opening all flight NetCDF Files
This file contains PRELIMINARY DATA that are NOT to be used for critical analysis.
NIDAS version: v1.2.3-115
NetCDF: Attribute not found
Processing Date & Time: 2024-11-11T22:04:13 +0000
Done reading ff01
Done reading ff02
Done reading ff03
Done reading ff04
Done reading ff05
Done reading ff06
Done reading rf01
Done reading rf02
Done reading rf03
Done reading rf04
Done reading rf05
Done reading rf06
Done reading rf07
Done reading rf08
Done reading rf09
Done reading rf10
Done reading tf01
Done reading tf02
Done reading tf03
Done reading tf04
Done reading tf05
Done reading tf06

Vertical Velocity Validations

Pitch Maneuvers

• ARISTO (2015) coefficients: 4.7523, 9.7908, 6.0781
• WECAN (2019) coefficients: 5.3236, 13.1197, 3.1606
• Test Coefficients: 5.4144, 13.3328, 3.6114

TF01

importlib.reload(caesar_utils)
flight = 'tf01'
mask = caesar_utils.mask_in_times(data[flight], *pitch_times[flight][0])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_time_series_pitch(data[flight], data[flight], mask, mask)

TF02 Leg 1

importlib.reload(caesar_utils)
flight = 'tf02'
mask = caesar_utils.mask_in_times(data[flight], *pitch_times[flight][0])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_time_series_pitch(data[flight], data[flight], mask, mask)

TF02 Leg 2

importlib.reload(caesar_utils)
flight = 'tf02'
mask = caesar_utils.mask_in_times(data[flight], *pitch_times[flight][1])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_time_series_pitch(data[flight], data[flight], mask, mask)

TF05 Leg 1

importlib.reload(caesar_utils)
flight = 'tf05'
mask = caesar_utils.mask_in_times(data[flight], *pitch_times[flight][0])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_time_series_pitch(data[flight], data[flight], mask, mask)

TF05 Leg 2

importlib.reload(caesar_utils)
flight = 'tf05'
mask = caesar_utils.mask_in_times(data[flight], *pitch_times[flight][1])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_time_series_pitch(data[flight], data[flight], mask, mask)

# compute metrics for each maneuver
importlib.reload(caesar_utils)
# first, mask data for each maneuver
maneuver_objs = []
for flight, times in pitch_times.items():
    for i, st_sp in enumerate(times):
        beg_time = st_sp[0]
        end_time = st_sp[1]
        title = f"Speed run {i+1} on {flight}"
        maneuv_mask = caesar_utils.mask_in_times(data[flight], beg_time=beg_time, end_time=end_time)
        fl_obj = caesar_utils.aoa_fit(data[flight], maneuv_mask, flight, leg=(i+1))
        maneuver_objs.append(fl_obj)

# iterate over all maneuver objects
for obj in maneuver_objs:
    aoa_bias = np.mean(obj.akrd - obj.aoa_ref)
    vspd_std = np.std(obj.vspd)
    w_std = np.std(obj.w)
    w_frac = w_std/vspd_std*100
    print(f"flight: {obj.flight}, aoa_bias: {aoa_bias:5.2f}, vspd_std: {vspd_std:4.2f}, w_std: {w_std:4.2f}, w_frac: {w_frac:5.2f}%")

flight: tf01, aoa_bias: -0.02, vspd_std: 5.90, w_std: 0.25, w_frac:  4.23%
flight: tf02, aoa_bias:  0.27, vspd_std: 5.30, w_std: 1.24, w_frac: 23.33%
flight: tf02, aoa_bias:  0.10, vspd_std: 4.20, w_std: 0.87, w_frac: 20.79%
flight: tf05, aoa_bias:  0.05, vspd_std: 4.01, w_std: 0.15, w_frac:  3.82%
flight: tf05, aoa_bias: -0.02, vspd_std: 4.41, w_std: 0.15, w_frac:  3.38%

W as a function of altitude

bdifr_blank_outs = {'rf05': [['06:53:00', '10:14:21']]}
adifr_blank_outs = {'rf02': [['17:54:50', '17:58:32']], 'rf07': [['15:56:17', '16:02:53'], ['16:10:50', '16:19:10']]}
adifr_icing = {'rf06': [['07:45:16', '07:50:12'], 
                        ['07:56:55', '08:05:06'], 
                        ['10:15:13', '10:21:58'], 
                        ['10:26:10', '11:02:31']], 
               'rf01': [['16:09:31', '16:36:16']], 
               'rf02': [['15:56:30', '16:02:12'], 
                        ['16:43:18', '18:16:00']], 
               'rf04': [['09:45:00', '09:48:40']], 
               'rf05': [['09:45:09', '10:16:58'], 
                        ['10:19:27', '10:35:08'], 
                        ['10:47:26', '11:22:03']], 
               'rf07': [['14:43:57', '15:26:56'], 
                        ['15:39:29', '15:56:17'], 
                        ['16:19:10', '17:35:00']], 
               'rf09': [['10:10:50', '10:52:10'], 
                        ['13:13:49', '14:23:30'], 
                        ['14:29:46', '14:45:54'], 
                        ['14:45:54', '16:02:00']], 
               'rf10': [['11:12:56', '11:20:58'], 
                        ['11:22:37', '11:27:59'], 
                        ['13:36:38', '13:40:27'], 
                        ['14:16:56', '14:30:50']], 
               'tf06': [['11:12:56', '11:20:58'], 
                        ['11:22:37', '11:27:59'], 
                        ['13:36:38', '13:40:27'], 
                        ['14:16:56', '14:30:50']]}


# concatenating into one dictionary
all_blank_outs = copy.deepcopy(adifr_icing)
all_blank_outs['rf05'].append(bdifr_blank_outs['rf05'][0])
all_blank_outs['rf02'].append(adifr_blank_outs['rf02'][0])
all_blank_outs['rf07'].append(adifr_blank_outs['rf07'][0])
all_blank_outs['rf07'].append(adifr_blank_outs['rf07'][1])
# converting hms dictionary to sfm dictionary
all_blank_outs_sfm = caesar_utils.hms_dict_to_sfm(all_blank_outs)

importlib.reload(caesar_utils)
# iterate through all flights, collect all flight objects
fl_objs = []
for flight, df in data.items():
    #mask = caesar_utils.mask_flying(df)
    mask = caesar_utils.mask_straight_and_level(df)
    no_icing_mask = mask
    if flight in all_blank_outs_sfm:
        for pair in all_blank_outs_sfm[flight]:
            no_icing_mask = np.logical_and(no_icing_mask, caesar_utils.mask_out_times(df, *pair))
    fl_objs.append(caesar_utils.aoa_fit(df, no_icing_mask, flight=flight))

# append all flight objects
all_data = copy.deepcopy(fl_objs[0])
for fl_obj in fl_objs[1:]:
    all_data = all_data.append(fl_obj)

caesar_utils.plot_all_z_vs_w(fl_objs, n_z_bins=60)

importlib.reload(caesar_utils)
caesar_utils.plot_z_vs_w(all_data, n_z_bins=60)

min z: 152.9105682373047, max z: 8247.2919921875, bin_depth: 134.90635681152344

Horizontal Wind Validations

Here, oscillations in $u$ and $v$ are compared against the lateral wind in the aircraft reference frame, estimated from the reference angle-of-sideslip.

The reference angle-of-sideslip, $\beta^{*}$, is derived from the first-order equations for $u$ and $v$:

\begin{equation} u = -V\sin(\Psi + \beta) + u_{p} \end{equation}

\begin{equation} v = -V\cos(\Psi + \beta) + v_{p} \end{equation}

Solving for $\beta$:

\begin{equation} \beta = -\Psi + \tan^{-1}\left(\frac{u_{p}-u}{v_{p}-v}\right) \end{equation}

Unfortunately, this equation for $\beta$ depends on $u$ and $v$, which depends on $\beta$. What is typically done is $u$ and $v$, previously estimated using previous calibration coefficients, are smoothed to remove any oscillations associated with the yaw maneuver. Here, this was done using a 5th Butterworth digital IIR filter with a cutoff period of 30 s to remove the yaw maneuvers with a period of $\approx15$ s. This filter was applied forward and backward to prevent phase shifting. With overbars indicating low-pass smoothing, the reference angle of sideslip is

\begin{equation} \beta^{*} = -\Psi + \tan^{-1}\left(\frac{u_{p}-\overline{u}}{v_{p}-\overline{v}}\right) \end{equation}

This reference angle of sideslip can then be used to estimate the lateral speed of the aircraft through the air which is being accounted for and removed from the horizontal wind components:

Lateral wind speed/TAS: \begin{equation} v_{y}^{*} = V\sin\beta^{*} \end{equation}

If variations in $u$, $v$ are small wrt $v_{y}^{*}$, then the calibration is considered successfull.

Yaw Maneuvers

• ARISTO (2016) coefficients: 1.5478, 12.3612
• WECAN (2019) coefficients: 0.85, 12.6582
• Test Coefficients: 1.6087, 13.2255

TF01

importlib.reload(caesar_utils)
all_curr_objs = {}
for flight, pairs in yaw_times.items():
    for i, pair in enumerate(pairs):
        maneuv_mask_prev = caesar_utils.mask_in_times(data[flight], *pair)
        curr_obj = caesar_utils.aos_fit(data[flight], maneuv_mask_prev, flight)
        key = f"{flight}_{i}"
        all_curr_objs[key] = curr_obj

importlib.reload(caesar_utils)
leg = 0
flight = 'tf01'
key = f"{flight}_{leg}"
mask = caesar_utils.mask_in_times(data[flight], *yaw_times[flight][leg])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_aos_validation(all_curr_objs[key], all_curr_objs[key])

TF02 Leg 1

importlib.reload(caesar_utils)
flight = 'tf02'
leg = 0
key = f"{flight}_{leg}"
mask = caesar_utils.mask_in_times(data[flight], *yaw_times[flight][leg])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_aos_validation(all_curr_objs[key], all_curr_objs[key])

TF02 Leg 2

importlib.reload(caesar_utils)
flight = 'tf02'
leg = 1
key = f"{flight}_{leg}"
mask = caesar_utils.mask_in_times(data[flight], *yaw_times[flight][leg])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_aos_validation(all_curr_objs[key], all_curr_objs[key])

TF05 Leg 1

importlib.reload(caesar_utils)
flight = 'tf05'
leg = 0
key = f"{flight}_{leg}"
mask = caesar_utils.mask_in_times(data[flight], *yaw_times[flight][leg])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_aos_validation(all_curr_objs[key], all_curr_objs[key])

TF05 Leg 2

importlib.reload(caesar_utils)
flight = 'tf05'
leg = 1
key = f"{flight}_{leg}"
mask = caesar_utils.mask_in_times(data[flight], *yaw_times[flight][leg])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_aos_validation(all_curr_objs[key], all_curr_objs[key])

TF05 Leg 3

importlib.reload(caesar_utils)
flight = 'tf05'
leg = 2
key = f"{flight}_{leg}"
mask = caesar_utils.mask_in_times(data[flight], *yaw_times[flight][leg])
caesar_utils.plot_track(data[flight], mask, title=flight)

# plot pitch maneuvers, see if aircraft motion is retained
importlib.reload(caesar_utils)
caesar_utils.plot_aos_validation(all_curr_objs[key], all_curr_objs[key])

# roughly quantify amount of lateral wind speed remaining in wind traces
importlib.reload(caesar_utils)
for key in all_curr_objs.keys():
    # used high-pass filtered u and v
    c_u_hp = all_curr_objs[key].u - all_curr_objs[key].u_sm
    c_v_hp = all_curr_objs[key].v - all_curr_objs[key].v_sm
    # for denominator, used the lateral speed from the reference angle-of-attack
    c_vy = all_curr_objs[key].lat_spd

    c_u_std = np.std(c_u_hp)
    c_v_std = np.std(c_v_hp)
    c_vy_std = np.std(c_vy)
    
    c_u_frac = c_u_std/c_vy_std*100
    c_v_frac = c_v_std/c_vy_std*100
    #print(f"leg: {key}, Prev u frac: {p_u_frac:6.3f}, v frac: {p_v_frac:6.3f}")
    print(f"leg: {key}, Curr u frac: {c_u_frac:6.3f}%, v frac: {c_v_frac:6.3f}%")

leg: tf01_0, Curr u frac:  4.223%, v frac:  3.168%
leg: tf02_0, Curr u frac: 10.826%, v frac: 11.480%
leg: tf02_1, Curr u frac: 11.493%, v frac: 15.768%
leg: tf05_0, Curr u frac:  2.214%, v frac:  2.622%
leg: tf05_1, Curr u frac:  2.957%, v frac:  2.033%
leg: tf05_2, Curr u frac:  4.537%, v frac:  3.676%

Yaw Maneuver Summary

• Except for TF02 (turbulent), the ratios of $\sigma_{u}/\sigma_{u_{y}}$ and $\sigma_{v}/\sigma_{v_{y}}$ are less than 5%
• Variations in the demeaned heading time series are not completely removed from the wind direction, U, V time series
• However, the oscillations that remain are less than $O(1)$ m/s, while the lateral wind speeds inferred from the reference angle-of-sideslip are $O(10)$ m/s